Install necessary packages: Set working directory
library(ggplot2)
library(readr)
library(sf)
Linking to GEOS 3.9.3, GDAL 3.5.2, PROJ 8.2.1; sf_use_s2() is TRUE
library(rnaturalearth)
The legacy packages maptools, rgdal, and rgeos, underpinning the sp package,
which was just loaded, will retire in October 2023.
Please refer to R-spatial evolution reports for details, especially
https://r-spatial.org/r/2023/05/15/evolution4.html.
It may be desirable to make the sf package available;
package maintainers should consider adding sf to Suggests:.
The sp package is now running under evolution status 2
(status 2 uses the sf package in place of rgdal)
Support for Spatial objects (`sp`) will be deprecated in {rnaturalearth} and will be removed in a future release of the package. Please use `sf` objects with {rnaturalearth}. For example: `ne_download(returnclass = 'sf')`
library(rnaturalearthdata)
Attaching package: 'rnaturalearthdata'
The following object is masked from 'package:rnaturalearth':
countries110
library(rgeos)
Loading required package: sp
rgeos version: 0.6-3, (SVN revision 696)
GEOS runtime version: 3.9.3-CAPI-1.14.3
Please note that rgeos will be retired during October 2023,
plan transition to sf or terra functions using GEOS at your earliest convenience.
See https://r-spatial.org/r/2023/05/15/evolution4.html for details.
GEOS using OverlayNG
Linking to sp version: 2.0-0
Polygon checking: TRUE
library(raster)
library(scales)
Attaching package: 'scales'
The following object is masked from 'package:readr':
col_factor
library(tmap)
setwd('C:/Users/cy_su/PycharmProjects/DSCI_605_Data_Visualizations/Module 6/M6_Lab5/')
Lab 1: Load the Data:
plot_locations <- read_csv("HARV_PlotLocations.csv")
Rows: 21 Columns: 16
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
chr (11): geodeticDa, utmZone, plotID, stateProvi, county, domainName, domai...
dbl (5): easting, northing, plotSize, elevation, plotdim_m
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
raster_image <- raster("HARV_chmCrop.tif")
Lab 1: Convert raster to a data frame for ggplot2:
raster_df <- as.data.frame(rasterToPoints(raster_image))
names(raster_df) <- c("easting", "northing", "value")
Lab 1: Plot the Raster Image: Provide a legend title and better color gradient for the raster. Add an outline to the points and makes them more visible. Add a caption and refines labels. Apply a minimalistic theme for a cleaner look. Adjust text sizes and styles for better readability and aesthetics.
ggplot() +
geom_raster(data = raster_df, aes(x = easting, y = northing, fill = value)) +
scale_fill_gradientn(colors = terrain.colors(10), name = "Elevation") +
geom_point(data = plot_locations, aes(x = easting, y = northing), color = "red", size = 4, shape = 21, fill = "white") +
labs(title = "Plot Locations on Elevation Map", x = "Easting (m)", y = "Northing (m)", caption = "Source: HARV Dataset") +
theme_minimal() +
theme(
plot.title = element_text(face = "bold", hjust = 0.5, size = 20),
axis.title = element_text(face = "bold", size = 14),
legend.title = element_text(face = "bold", size = 12),
legend.text = element_text(size = 12),
plot.caption = element_text(hjust = 0, size = 10)
) +
coord_cartesian(expand = FALSE)
Lab 2: Load the Data:
watersheds <- st_read("Watersheds_HUC08_2009/WATERSHEDS_HUC08_2009_USDA_IN.shp")
Reading layer `WATERSHEDS_HUC08_2009_USDA_IN' from data source
`C:\Users\cy_su\PycharmProjects\DSCI_605_Data_Visualizations\Module 6\M6_Lab5\Watersheds_HUC08_2009\WATERSHEDS_HUC08_2009_USDA_IN.shp'
using driver `ESRI Shapefile'
Simple feature collection with 39 features and 12 fields
Geometry type: POLYGON
Dimension: XY
Bounding box: xmin: 374222.5 ymin: 4157935 xmax: 794277.8 ymax: 4691819
Projected CRS: NAD83 / UTM zone 16N
states <- st_read("tl_2019_us_state/tl_2019_us_state.shp")
Reading layer `tl_2019_us_state' from data source
`C:\Users\cy_su\PycharmProjects\DSCI_605_Data_Visualizations\Module 6\M6_Lab5\tl_2019_us_state\tl_2019_us_state.shp'
using driver `ESRI Shapefile'
Simple feature collection with 56 features and 14 fields
Geometry type: MULTIPOLYGON
Dimension: XY
Bounding box: xmin: -179.2311 ymin: -14.60181 xmax: 179.8597 ymax: 71.43979
Geodetic CRS: NAD83
Set the theme:
tmap_mode("view")
tmap mode set to interactive viewing
Plot the Data: Set a color palette. Change the color of states to a simple grey to make the watersheds more visible For the watersheds adds a white border and sets up a pop-up showing the square miles of the area when an area is clicked on
palette <- RColorBrewer::brewer.pal(9, "YlGnBu")
tm_shape(states) +
tm_polygons(col = "grey70", border.col = "grey40") +
tm_shape(watersheds) +
tm_polygons(col = "HUC_8", palette = palette, title = "Watershed (HUC_8)",
border.col = "white", border.alpha = .5, id = "HUC_8",
popup.vars = c("Area (sq miles)" = "SQ_MILES")) +
tm_compass(type = "arrow", position = c("left","top"), size = 1.5) +
tm_scale_bar(breaks = c(0, 50, 100), position = c("right","bottom"), size = 1) +
tm_layout(main.title = "Watershed Map of Indiana (HUC_8)",
main.title.position = "center",
legend.position = c("left","bottom"),
legend.text.size = 0.7,
legend.title.size = 0.9,
legend.bg.color = "white",
legend.bg.alpha = 0.5,
legend.frame = TRUE,
frame = FALSE)
Warning: The argument size of tm_scale_bar is deprecated. It has been renamed to
text.size
Compass not supported in view mode.
legend.postion is used for plot mode. Use view.legend.position in tm_view to set the legend position in view mode.
Warning: Number of levels of the variable "HUC_8" is 39, which is
larger than max.categories (which is 30), so levels are combined. Set
tmap_options(max.categories = 39) in the layer function to show all levels.
Warning: In view mode, scale bar breaks are ignored.